
! This contains generic routines
! Random number routines here are scheduled for withdrawal. Need to replace with new ones
    module routines_generic
      use precisions

    contains
    
    
      subroutine linearinterp1(xgrid, ygrid, length, xval, yval, extrapbot, extraptop)
! Subroutine for linear interpolation
! Arguments are 1. (input) xvlalues
! 2. (input) yvalues
! 3. (input) length of xgrid and ygrid,
! 4. (input) value of x to interpolate at
! 5. (output) interpolated value
! 6. (input) indicator for whether to extrapolate below bottom point of grid
! 7. (input indicator for whether to extrapolate above top point of grid

! If the choose not to extrapolate interpolated value is set equal to either the top or the bottom
! point of the grid
! It is fully generic

        implicit none
! Arguments
        integer, intent (in) :: length
        integer, intent (in) :: extrapbot
        integer, intent (in) :: extraptop
        real (kind=rk), intent (in) :: xgrid(length)
        real (kind=rk), intent (in) :: ygrid(length)
        real (kind=rk), intent (in) :: xval
        real (kind=rk), intent (out) :: yval

! For programme
        integer :: xlowerloc
        integer scratch_xnearestloc ! This areguments required by findingridr but I'm not using them here - hence
! scratch

        intrinsic size

! Checking inputs - have commented this out for speed
     !   if ((size(xgrid).ne.length) .or. (size(ygrid).ne.length)) then
     !     write (*, *) 'An argument in linearinterp1 is not of the specified length'
     !     stop
     !   end if

     !   if ((extrapbot.ne.0 .and. extrapbot.ne.1) .or. (extraptop.ne.0 .and. extraptop.ne.1)) then
     !     write (*, *) 'An argument in linearterp1 that needs to be either 0 or 1 is not so'
     !     stop
     !   end if

        if (abs(xgrid(length)-xgrid(1)).lt.0.1d-15) then ! If the top of the grid is equal to the bottom
          yval = ygrid(1) ! (which would include the scenario that all the entries are the same) then I
! assign the bottom value of y to yval
        else

          call findingridr(xgrid, length, xval, scratch_xnearestloc, xlowerloc)
          call linearinterpfromlocations(xgrid, ygrid, xlowerloc, xval, length, extrapbot, extraptop, yval)
        end if

      end subroutine linearinterp1

! ---------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------!


! Extrapolation is by default here. There is no option to use the bottom value (as I have
! in the one-dimensional linear interpolation).

      subroutine linearinterp2_withextrap(x1grid, x2grid, dim1, dim2, x1val, x2val, yval, extrapbot, extraptop, ygrid)
        implicit none
! Arguments
        integer, intent (in) :: dim1
        integer, intent (in) :: dim2
        real (kind=rk), intent (in) :: x1grid(dim1)
        real (kind=rk), intent (in) :: x2grid(dim2)
        real (kind=rk), intent (in) :: ygrid(dim1, dim2)
        real (kind=rk), intent (in) :: x1val
        real (kind=rk), intent (in) :: x2val
        real (kind=rk), intent (out) :: yval
        integer, intent (in) :: extrapbot
        integer, intent (in) :: extraptop

! For programme
        integer :: x1lowerloc
        integer :: x2lowerloc
        integer :: scratch_nearestloc
        integer :: newxgridlowerloc
        real (kind=rk) :: newygrid(2)
        real (kind=rk) :: newxgrid(2)

        if ((size(x1grid).ne.dim1) .or. (size(x2grid).ne.dim2)) then
          write (*, *) 'An argument in linearinterp2 is not of the specified length'
          stop
        end if

! Find the square
        call findingridr(x1grid, dim1, x1val, scratch_nearestloc, x1lowerloc)
        call findingridr(x2grid, dim2, x2val, scratch_nearestloc, x2lowerloc)

! Find if any of the values are off the grid
        if (x1lowerloc.eq.0) then
          x1lowerloc = 1
        end if

        if (x2lowerloc.eq.0) then
          x2lowerloc = 1
        end if

        if (x1lowerloc.eq.dim1) then
          x1lowerloc = dim1 - 1
        end if

        if (x2lowerloc.eq.dim2) then
          x2lowerloc = dim2 - 1
        end if

! Now do linearinterpolation along the first dimension
        call linearinterpfromlocations(x1grid, ygrid(:,x2lowerloc), x1lowerloc, x1val, dim1, extrapbot, extraptop, &
          newygrid(1))
        call linearinterpfromlocations(x1grid, ygrid(:,(x2lowerloc+1)), x1lowerloc, x1val, dim1, extrapbot, extraptop, &
          newygrid(2))

! Now generate a new grid from from the x2 bounds (a new ygrid has already bee generated)
        newxgrid(1) = x2grid(x2lowerloc)
        newxgrid(2) = x2grid(x2lowerloc+1)

        if (x2lowerloc.eq.0) then
          newxgridlowerloc = 0
        else if (x2lowerloc.eq.dim2) then
          newxgridlowerloc = 2
        else
          newxgridlowerloc = 1
        end if
        call linearinterpfromlocations(newxgrid, newygrid, newxgridlowerloc, x2val, 2, extrapbot, extraptop, yval)

      end subroutine linearinterp2_withextrap


! ---------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------!
! Extrapolation is by default here. There is no option to use the bottom value (as I have
! in the one-dimensional linear interpolation).

      subroutine linearinterp2(x1grid, x2grid, dim1, dim2, x1val, x2val, yval, ygrid)
!        use types
        implicit none
! Arguments
        integer, intent (in) :: dim1
        integer, intent (in) :: dim2
        real (kind=rk), intent (in) :: x1grid(dim1)
        real (kind=rk), intent (in) :: x2grid(dim2)
        real (kind=rk), intent (in) :: x1val
        real (kind=rk), intent (in) :: x2val
        real (kind=rk), intent (out) :: yval
        real (kind=rk), intent (in) :: ygrid(dim1, dim2)
        
! For programme
        integer :: x1lowerloc
        integer :: x2lowerloc
        integer :: scratch_nearestloc
        integer :: newxgridlowerloc
        real (kind=rk) :: newygrid(2)
        real (kind=rk) :: newxgrid(2)
        
        !Have commented out for speed
        if ((size(x1grid).ne.dim1) .or. (size(x2grid).ne.dim2)) then
          write (*, *) 'An argument in linearinterp2 is not of the specified length'
          stop
        end if

        call findingridr(x1grid, dim1, x1val, scratch_nearestloc, x1lowerloc)

        if (x1lowerloc.eq.0) then
          x1lowerloc = 1
        end if
        if (x1lowerloc.eq.dim1) then
          x1lowerloc = dim1 - 1
        end if

! Find the square
          call findingridr(x2grid, dim2, x2val, scratch_nearestloc, x2lowerloc)

! Find if any of the values are off the grid
        
          if (x2lowerloc.eq.0) then
            x2lowerloc = 1
          end if
          if (x2lowerloc.eq.dim2) then
            x2lowerloc = dim2 - 1
          end if



         
         !write(*,*) x1lowerloc
         !write(*,*) x2lowerloc
         !write(*,*) dim1 
         !write(*,*) ygrid(:,x2lowerloc)
         !write(*,*) x1grid
         
! Now do linearinterpolation along the first dimension
!        if (check1) then
!        write(*,*) x1lowerloc, x2lowerloc
!        write(*,*) dim1, dim2
!        write(*,*) x1val 
!        write(*,*) ygrid
!        stop
!        end if

          call linearinterpfromlocations(x1grid, ygrid(:,x2lowerloc), x1lowerloc, x1val, dim1, 1, 1, newygrid(1))
          call linearinterpfromlocations(x1grid, ygrid(:,(x2lowerloc+1)), x1lowerloc, x1val, dim1, 1, 1, newygrid(2))


! Now generate a new grid from from the x2 bounds (a new ygrid has already bee generated)
        newxgrid(1) = x2grid(x2lowerloc)
        newxgrid(2) = x2grid(x2lowerloc+1)



!Think this is redundant
        if (x2lowerloc.eq.0) then
          newxgridlowerloc = 0
        else if (x2lowerloc.eq.dim2) then
          newxgridlowerloc = 2
        else
          newxgridlowerloc = 1
        end if
           
           



        call linearinterpfromlocations(newxgrid, newygrid, newxgridlowerloc, x2val, 2, 1, 1, yval)

      end subroutine linearinterp2


! ---------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------!

! Subroutine for interpolating in one dimension given the points on the grid that bound the value for
! evaluation
      subroutine linearinterpfromlocations(xgrid, ygrid, xlowerloc, xval, length, extrapbot, extraptop, yval)
        implicit none
        integer, intent (in) :: length
        real (kind=rk), intent (in) :: xgrid(length)
        real (kind=rk), intent (in) :: ygrid(length)
        integer, intent (in) :: xlowerloc
        real (kind=rk), intent (in) :: xval
        integer, intent (in) :: extrapbot
        integer, intent (in) :: extraptop
        real (kind=rk), intent (out) :: yval


! For programme
        real (kind=rk) xgridmin
        real (kind=rk) xgridmax
        real (kind=rk) xgriddiffinv
        real (kind=rk) weightmin
        real (kind=rk) weightmax
        real (kind=rk) extrapslope
        integer :: xupperloc

! First, and easiest, get the interpolated value if xval is interior to the grid
        if (xlowerloc.ne.0 .and. xlowerloc.ne.length) then
          xupperloc = xlowerloc + 1
          xgridmin = xgrid(xlowerloc)
          xgridmax = xgrid(xupperloc)
          xgriddiffinv = 1/(xgridmax-xgridmin)

          weightmax = (xval-xgridmin)*xgriddiffinv
          weightmin = (xgridmax-xval)*xgriddiffinv

! Calculating the interpolated function values
          yval = weightmin*ygrid(xlowerloc) + weightmax*ygrid(xupperloc)
        else if (xlowerloc.eq.0) then
          if (extrapbot.eq.0) then
            yval = ygrid(1)
          else
            extrapslope = (ygrid(2)-ygrid(1))/(xgrid(2)-xgrid(1))
            yval = ygrid(1) + extrapslope*(xval-xgrid(1))
          end if

        else ! a is hu-eg so we extrapolate
          if (extraptop.eq.0) then
            yval = ygrid(length)
          else
            extrapslope = (ygrid(length)-ygrid(length-1))/(xgrid(length)-xgrid(length-1))
            yval = ygrid(length) + extrapslope*(xval-xgrid(length))
          end if
        end if


      end subroutine linearinterpfromlocations

! ---------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------!

! ---------------------------------------------------------------------------------------!
!------------------------------------------------------------------------!
! Fully generic - need to write explanation
      
      
      subroutine gengrid_generic(lowerbounds, upperbounds, numyears, gridpoints, growthgrid, assgrid)
        implicit none

! Arguments
        integer, intent (in) :: numyears
        integer, intent (in) :: gridpoints
        real (kind=rk), intent (in) :: lowerbounds(numyears)
        real (kind=rk), intent (in) :: upperbounds(numyears)
        real (kind=rk), intent (in) :: growthgrid

        real (kind=rk), intent (out) :: assgrid(gridpoints, (numyears+1))

! For program
        integer ixassnode, ixtoday
       

! Getting the asset grid
        do ixtoday = numyears, 1, -1
          call getnodes(assgrid(:, ixtoday), lowerbounds(ixtoday), upperbounds(ixtoday), gridpoints, growthgrid)
        end do ! time loop

      end subroutine gengrid_generic
! ---------------------------------------------------------------------------------------!
! Fully generic - need to write explanation

      subroutine gengrid_extra(lowerbounds, breaks1, breaks2, upperbounds, numyears, gridpoints1, gridpoints2, &
          gridpoints3, growthgrid1, growthgrid2, growthgrid3, gridpoints, assgrid)
        implicit none

! Arguments
        integer, intent (in) :: numyears
        integer, intent (in) :: gridpoints1
        integer, intent (in) :: gridpoints2
        integer, intent (in) :: gridpoints3
        integer, intent (in) :: gridpoints
        real (kind=rk), intent (in) :: lowerbounds(numyears)
        real (kind=rk), intent (in) :: breaks1(numyears)
        real (kind=rk), intent (in) :: breaks2(numyears)
        real (kind=rk), intent (in) :: upperbounds(numyears)
        real (kind=rk), intent (in) :: growthgrid1
        real (kind=rk), intent (in) :: growthgrid2
        real (kind=rk), intent (in) :: growthgrid3

        real (kind=rk), intent (out) :: assgrid(gridpoints, numyears)

! For program
        integer ixassnode, ixtoday, index
        real (kind=rk) restofgrid1(gridpoints1)
        real (kind=rk) restofgrid2(gridpoints2+1)
        real (kind=rk) restofgrid3(gridpoints3+1)

! Check gridpoints
        if ((gridpoints1+gridpoints2+gridpoints3).ne.gridpoints) then
          write (*, *) 'Error in gengrid_extra; gridpoints1 + gridpoints2 + gridpoints3 is not equal gridpoints'
          stop
        end if

! Getting the asset grid
        do ixtoday = 1, numyears, 1

! I start the grid from the constraint plus the mininmum shock
          !assgrid(1, ixtoday) = lowerbounds(ixtoday)
          !assgrid(2, ixtoday) = assgrid(1, ixtoday) + 20.0_rk

          call getnodes(restofgrid1, lowerbounds(ixtoday), breaks1(ixtoday), gridpoints1, growthgrid1)
          do ixassnode = 1, gridpoints1 - 1, 1 ! Don't want to add in break1 - leave that till next one
            assgrid(ixassnode, ixtoday) = restofgrid1(ixassnode)
          end do

          call getnodes(restofgrid2, breaks1(ixtoday), breaks2(ixtoday), gridpoints2+1, growthgrid2)
          index = 1
          do ixassnode = (gridpoints1), (gridpoints1+gridpoints2-1), 1
            assgrid(ixassnode, ixtoday) = restofgrid2(index)
            index = index + 1
          end do

          call getnodes(restofgrid3, breaks2(ixtoday), upperbounds(ixtoday), gridpoints3+1, growthgrid3)
          index = 1
          do ixassnode = (gridpoints1+gridpoints2), gridpoints, 1
            assgrid(ixassnode, ixtoday) = restofgrid3(index)
            index = index + 1
          end do

        end do ! time loop

      end subroutine gengrid_extra
! ---------------------------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------------------------!

! This function returns a grid where:
! 1. The grid starts at firstpoint (input)
! 2. The grid ends at lastpoint (input)
! 3. The number of nodes is numberofnodes (input)
! 4. The growth rate of the step size is growthrate (input)

! This function is fully generic

      subroutine getnodes(nodes, firstpoint, lastpoint, numberofnodes, growthrate, logspace)
        implicit none

! Arguments
        real (kind=rk), intent (in) :: firstpoint
        real (kind=rk), intent (in) :: lastpoint
        real (kind=rk), intent (in) :: growthrate
        integer, intent (in) :: numberofnodes
        real (kind=rk), intent (out) :: nodes(numberofnodes)
        logical, optional, intent(in) ::  logspace 

! For programme
        real (kind=rk) step, span, logStep, lognodes(numberofnodes)
        real (kind=rk) denom
        integer ix
        !intrinsic :: log
  
        
! Check inputs
        if (firstpoint.ge.lastpoint) then
          write (*, *) firstpoint, lastpoint
          write (*, *) 'In getnodes, firstpoint must be strictly less than lastpoint'
          write (*, *) 'The +10._rk is clumsy in gengrid_Extra'
          stop
        end if

! First I need, conditional, on the span of the space, the number of gridpoints, and the growth
! rate of the space between successive gridpoints, to work out the size of the first step
        denom = 0
        do ix = 1, numberofnodes - 1, 1
          denom = denom + (1+growthrate)**(ix-1)
        end do

        step = (lastpoint-firstpoint)/denom
        
        if (present(logspace)) then
            if (.not.(logspace)) then
                write(*,*) 'only include logspace in getnodes if you intend it to be true'
                stop
            end if
            if (numberofnodes.eq.1) then
                write(*,*) 'cant use one node here as friends dont let friends divide by zero'
                stop
            end if
            
            span = lastpoint-firstpoint
            logStep = (log(lastpoint+1) - log(firstpoint+1))/(numberofnodes-1)
            lognodes(1) = log(firstpoint+1)
             do ix = 2, numberofnodes, 1
                lognodes(ix) = lognodes(ix-1) + logStep
             end do
             nodes = exp(lognodes) -( firstpoint+1)
            else
                nodes(1) = firstpoint
                do ix = 2, numberofnodes, 1
                  nodes(ix) = nodes(ix-1) + step
                  step = step*(1+growthrate)
                end do
            end if
            

      end subroutine getnodes
! ---------------------------------------------------------------------------------------------------------!
! Function to get numerical first derivative. Generic other than use of linearinterp1 and I've hardcoded
! last two arguments in linearinterp1
      real (kind=rk) function deriv(x, fx, x0, length)
        implicit none

! Arguments
        integer, intent (in) :: length
        real (kind=rk), intent (in) :: x(length)
        real (kind=rk), intent (in) :: fx(length)
        real (kind=rk), intent (in) :: x0

! For progamme
        real (kind=rk) :: forderiv1, forderiv2

        call linearinterp1_fast(x, fx, length, x0+0.0000001_rk, forderiv1, 0, 1)
        call linearinterp1_fast(x, fx, length, x0-0.0000001_rk, forderiv2, 0, 1)
        deriv = (forderiv1-forderiv2)/(0.0000002_rk)

      end function deriv
! ---------------------------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------------------------!
! Function to get the numerical derivative of a function when that function is defined on a 2 dimensional
! grid (with two arguments). Generic except for the reliance on linearinterp2
! NB: The derivative is with respect to the variable in x1grid at x1val
      real (kind=rk) function deriv_2d(x1grid, x2grid, ygrid, dim1, dim2, x1val, x2val)

        implicit none

! Arguments
        integer, intent (in) :: dim1
        integer, intent (in) :: dim2
        real (kind=rk), intent (in) :: x1grid(dim1)
        real (kind=rk), intent (in) :: x2grid(dim2)
        real (kind=rk), intent (in) :: ygrid(dim1, dim2)
        real (kind=rk), intent (in) :: x1val
        real (kind=rk), intent (in) :: x2val


! For progamme
        real (kind=rk) :: forderiv1, forderiv2

        !call linearinterp2(x1grid, x2grid, dim1, dim2, (x1val+0.0000001_rk), x2val, forderiv1, ygrid)
        !call linearinterp2(x1grid, x2grid, dim1, dim2, (x1val-0.0000001_rk), x2val, forderiv2, ygrid)
        call linearinterp2_fast(x1grid, x2grid, dim1, dim2, (x1val+0.0000001_rk), x2val, forderiv1, ygrid)
        call linearinterp2_fast(x1grid, x2grid, dim1, dim2, (x1val-0.0000001_rk), x2val, forderiv2, ygrid)
        deriv_2d = (forderiv1-forderiv2)/(0.0000002_rk)

      end function deriv_2d

! ---------------------------------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------------------------------!
! A subroutine for finding  the nodes and the values on a real grid that bound a value off the grid
! It is fully generic
      subroutine findingridr(xgrid, length, xval, nearestloc, lowerloc)
        implicit none
! Arguments
        integer, intent (in) :: length
        real (kind=rk), intent (in) :: xgrid(length)
        real (kind=rk), intent (in) :: xval

        integer, intent (out) :: nearestloc
        integer, intent (out) :: lowerloc

! For programme
        integer :: jl
        integer :: jm
        integer :: ju
        integer :: upperloc
        real (kind=rk) :: lowerval
        real (kind=rk) :: upperval
        intrinsic size

! Checking inputs
        if ((size(xgrid).ne.length)) then
          write (*, *) 'An argument in findingridr is not of the specified length'
          stop
        end if

! This block of code finds the two points on the xgrid between which the x value lies
        jl = 0
        ju = length + 1
        do while (ju-jl.gt.1)
          jm = (ju+jl)/2
          if (xval.gt.xgrid(jm)) then
            jl = jm
          else
            ju = jm
          end if
        end do

        lowerloc = jl
        upperloc = jl + 1

        if (lowerloc.eq.0) then
          lowerval = xgrid(1)
          upperval = xgrid(1)
        else if (lowerloc.eq.length) then
          lowerval = xgrid(length)
          upperval = xgrid(length)
        else
          lowerval = xgrid(lowerloc)
          upperval = xgrid(upperloc)
        end if

        if ((((upperval-xval).lt.(xval-lowerval)) .or. (lowerloc.eq.0)) .and. (lowerloc.ne.length)) then
          nearestloc = upperloc
        else
          nearestloc = lowerloc
        end if

        if ((nearestloc.lt.0) .or. (nearestloc.gt.length)) then
          write (*, *) 'error in findingr - nearestloc is outside the grid'
          stop
        end if
        
        if (nearestloc.eq.0) then
                write(*,*)
        end if
        

      end subroutine findingridr


      function transpose1dp(vector, length)
        implicit none
        integer, intent (in) :: length
        real (kind=rk), intent (in) :: vector(length)

        real (kind=rk) :: transpose1dp(1, length)

        integer :: ix

        do ix = 1, length, 1
          transpose1dp(1, ix) = vector(ix)
        end do

      end function transpose1dp
! ---------------------------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------------------------!
      function transpose1int(vector, length)
        implicit none
        integer, intent (in) :: length
        integer, intent (in) :: vector(length)

        integer :: transpose1int(1, length)

        integer :: ix

        do ix = 1, length, 1
          transpose1int(1, ix) = vector(ix)
        end do

      end function transpose1int
! ---------------------------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------------------------!


! ---------------------------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------------------------!
! A routine that takes a vector representing a PDF and returns a vector representing a CDF
      subroutine pdftocdf(length, pdf, cdf)
        implicit none
        integer, intent (in) :: length
        real (kind=rk), intent (in) :: pdf(length)
        real (kind=rk), intent (out) :: cdf(length)

! For programme
        integer :: ixlen
        real (kind=rk) :: sumforcdf

        if ((sum(pdf).lt.0.999) .or. (sum(pdf).gt.1.001)) then
          write (*, *) 'error in pdftocdf: probabilities in pdf do not sum to 1'
          write (*, *) 'sum is', sum(pdf)
          stop
        end if

! Get CDF
        sumforcdf = 0
        do ixlen = 1, length
          sumforcdf = sumforcdf + pdf(ixlen)
          cdf(ixlen) = sumforcdf
        end do


        if ((cdf(length).lt.1.0_rk).and.(cdf(length).gt.0.999)) then
            cdf(length) = 1.0_rk
        end if
        
      end subroutine pdftocdf
!---------------------------------------------------------------------------------------------------------!
!---------------------------------------------------------------------------------------------------------!
!Get uniform draws using the NAG libraries
      subroutine getdraws(length1, length2, totalrandom, seed, draws, normal, mu, var)
        implicit none

        integer, intent (in) :: length1, length2
        integer, intent (in) :: seed
        integer, intent (in) :: totalrandom
        logical, intent (in), optional :: normal
        real (kind=rk), intent (in), optional :: mu
        real (kind=rk), intent (in), optional :: var
        real (kind=rk), intent (out) :: draws(length1, length2)

!For programme
        integer :: genid
        integer :: subid !Says use the NAG basic generator
        integer :: lseed
        integer, parameter :: lstate = 633
        real (kind=rk) :: drawsvec(totalrandom)
        integer :: state(lstate)
        integer :: ifail
        external :: g05kff, g05saf

!Check that if we have requested normal variables, mu and var are provided
        if (present(normal)) then
          if ((normal) .and. (.not. (present(mu) .and. present(var)))) then
            write (*, *) 'error in getdraws. If normal is .true.,  mean and variance must be supplied'
            stop
          end if
        end if

!Check length1 * length2 = totalrandom
        if ((length1*length2).ne.totalrandom) then
          write (*, *) 'error in getunifdraws: length1 * length2 = totalrandom'
          stop
        end if

        genid = 1 !Says use the NAG basic generator; do not change - will not work otherwise; need to write new subroutine if want to use more genera
 !
!
!
!l
        subid = 1 !This is irrelavant when we're using the NAG basic generator (i.e. genid = 1)
        lseed = 1 !size of seed	
        ifail = 0

! Initialize the generator to a repeatable sequence
        call g05kff(genid, subid, seed, lseed, state, lstate, ifail)

        ifail = 0

        if (present(normal) .and. (normal)) then
          call g05skf(totalrandom, mu, var, state, drawsvec, ifail)
        else
          call g05saf(totalrandom, state, drawsvec, ifail)
        end if

        draws = reshape(drawsvec, (/length1,length2/) )


      end subroutine getdraws


! ---------------------------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------------------------!
!Function that gets the squared distance between two points
      function dist2(x, y)
        implicit none
        real (kind=rk), intent (in) :: x, y
        real (kind=rk) :: dist2

        dist2 = (x-y)*(x-y)

      end function dist2
!---------------------------------------------------------------------------------------------------------!
!---------------------------------------------------------------------------------------------------------!
! This is based on code given in Miranda and Fackler:
! "Applied Computaiontal Economics and Finance"

      real (kind=rk) function golden_generic(a, b, xmax, func)
        implicit none

! Generic stuff
! Arguments
        real (kind=rk), intent (in) :: a
        real (kind=rk), intent (in) :: b
        real (kind=rk) :: func
        real (kind=rk), intent (out) :: xmax

! For program
        real (kind=rk) :: x1
        real (kind=rk) :: x2
        real (kind=rk) :: f1
        real (kind=rk) :: f2
        real (kind=rk) :: diff
        real (kind=rk) :: goldenalpha1
        real (kind=rk) :: goldenalpha2
        real (kind=rk) :: tol
        parameter (goldenalpha2=0.611003399_rk)
        parameter (goldenalpha1=1.0_rk-goldenalpha2)
        parameter (tol=1.d-2)

        x1 = a + goldenalpha1*(b-a)
        x2 = a + goldenalpha2*(b-a)

        f1 = func(x1)
        f2 = func(x2)
        diff = goldenalpha1*goldenalpha2*(b-a)
        do while (diff.gt.tol)
          diff = diff*goldenalpha2
          if (f2.lt.f1) then
            x2 = x1
            x1 = x1 - diff

            f2 = f1
            f1 = func(x1)
          else
            x1 = x2
            x2 = x2 + diff

            f1 = f2
            f2 = func(x2)
          end if
        end do

        if (f2.gt.f1) then
          xmax = x2
          golden_generic = f2
        else
          xmax = x1
          golden_generic = f1
        end if

      end function golden_generic

!---------------------------------------------------------------------------------------------------------!
!---------------------------------------------------------------------------------------------------------!
!A subroutine for trilinear interpolation
!A nice method here (that I'm not using) method at: http://paulbourke.net/miscellaneous/interpolation/
!This has been tested with the subroutine test below
      subroutine linearinterp3(x1grid, x2grid, x3grid, dim1, dim2, dim3, x1val, x2val, x3val, yval, ygrid)
!        use types
        implicit none

! Arguments
        integer, intent (in) :: dim1
        integer, intent (in) :: dim2
        integer, intent (in) :: dim3
        real (kind=rk), intent (in) :: x1grid(dim1)
        real (kind=rk), intent (in) :: x2grid(dim2)
        real (kind=rk), intent (in) :: x3grid(dim3)
        real (kind=rk), intent (in) :: ygrid(dim1, dim2, dim3)
        real (kind=rk), intent (in) :: x1val
        real (kind=rk), intent (in) :: x2val
        real (kind=rk), intent (in) :: x3val
        real (kind=rk), intent (out) :: yval

! For programme
        integer :: x3lowerloc
        integer :: scratch_nearestloc
        integer :: newxgridlowerloc
        real (kind=rk) :: newygrid(2)
        real (kind=rk) :: newxgrid(2)


        if ((size(x1grid).ne.dim1) .or. (size(x2grid).ne.dim2) .or. (size(x3grid).ne.dim3)) then
          write (*, *) 'An argument in linearinterp3 is not of the specified length'
          stop
        end if

! Locate x3val on the grid
          call findingridr(x3grid, dim3, x3val, scratch_nearestloc, x3lowerloc)

! Find if any of the values are off the grid       
          if (x3lowerloc.eq.0) then
            x3lowerloc = 1
          end if

          if (x3lowerloc.eq.dim3) then
            x3lowerloc = dim3 - 1
          end if




!Do two two bilinear interpolations over x1grid and x2grid keeping  x3 at x3lowerloc and x3lowerloc + 1 
          call linearinterp2(x1grid, x2grid, dim1, dim2, x1val, x2val, newygrid(1), ygrid(:,:,x3lowerloc)) !fast
          call linearinterp2(x1grid, x2grid, dim1, dim2, x1val, x2val, newygrid(2), ygrid(:,:,(x3lowerloc+1)))

          ! call linearinterp2_fast(x1grid, x2grid, dim1, dim2, x1val, x2val, newygrid(1), ygrid(:,:,x3lowerloc))
          !call linearinterp2_fast(x1grid, x2grid, dim1, dim2, x1val, x2val, newygrid(2), ygrid(:,:,(x3lowerloc+1)))



! Now to a linear interpolation on these found vlues
! First generate a new grid from from the x2 bounds (a new ygrid has already bee generated)
        newxgrid(1) = x3grid(x3lowerloc)
        newxgrid(2) = x3grid(x3lowerloc+1)

!Think this is redundant
        if (x3lowerloc.eq.0) then
          newxgridlowerloc = 0
        else if (x3lowerloc.eq.dim3) then
          newxgridlowerloc = 2
        else
          newxgridlowerloc = 1
        end if

        call linearinterpfromlocations(newxgrid, newygrid, newxgridlowerloc, x3val, 2, 1, 1, yval)


        !if (yval.gt.100) then
        !    write(*,*) 't'
        !end if        

      end subroutine linearinterp3
!---------------------------------------------------------------------------------------------------------!
!---------------------------------------------------------------------------------------------------------!
!A subroutine for testing the routine above
      subroutine testlinearinterp3
        implicit none

!Number of test values
        integer, parameter :: numtests = 500

!How far outside the grid to check for extrapolation		(=1 for only check inside the grid)
        real (kind=rk), parameter :: outside = 1.0

!Dimensions
        integer, parameter :: x1dim = 50
        integer, parameter :: x2dim = 51
        integer, parameter :: x3dim = 52

!Bounds on the grids
        real (kind=rk), parameter :: low1 = 0.0, upp1 = 100.0
        real (kind=rk), parameter :: low2 = -100.0, upp2 = 300.0
        real (kind=rk), parameter :: low3 = -0.5, upp3 = 0.5

!Growth rate of the distance between gridpoints
        real (kind=rk), parameter :: growth1 = 0.05
        real (kind=rk), parameter :: growth2 = 0.00
        real (kind=rk), parameter :: growth3 = 0.00

!Grids
        real (kind=rk) :: x1grid(x1dim)
        real (kind=rk) :: x2grid(x2dim)
        real (kind=rk) :: x3grid(x3dim)
        real (kind=rk) :: ygrid(x1dim, x2dim, x3dim)

!Other stuff
        real (kind=rk), dimension (numtests, 1) :: x1val, x2val, x3val
        real (kind=rk), dimension (numtests) :: interpval, actval, absinterperror, relinterperror
        real (kind=rk) :: maxabserror, maxrelerror

!Indexes
        integer :: ix1, ix2, ix3, ixtest

!Fill in the x grids           
        call getnodes(x1grid, low1, upp1, x1dim, growth1)
        call getnodes(x2grid, low2, upp2, x2dim, growth2)
        call getnodes(x3grid, low3, upp3, x3dim, growth3)

!Fill in the y grids
        do ix3 = 1, x3dim, 1
          do ix2 = 1, x2dim, 1
            do ix1 = 1, x1dim, 1
              ygrid(ix1, ix2, ix3) = funcfortestinterp(x1grid(ix1), x2grid(ix2), x3grid(ix3))
            end do
          end do
        end do

!Get random numbers
        call getdraws(numtests, 1, numtests, 156898, x1val)
        x1val = (x1val*(upp1-low1)+low1)*outside

        call getdraws(numtests, 1, numtests, 256898, x2val)
        x2val = (x2val*(upp2-low2)+low2)*outside

        call getdraws(numtests, 1, numtests, 356898, x3val)
        x3val = (x3val*(upp3-low3)+low3)*outside

!Now do a test: get the exact value and the interpolated value for some points off the grid
        do ixtest = 1, numtests, 1
          call linearinterp3(x1grid, x2grid, x3grid, x1dim, x2dim, x3dim, x1val(ixtest,1), x2val(ixtest,1), &
            x3val(ixtest,1), interpval(ixtest), ygrid)
          actval(ixtest) = funcfortestinterp(x1val(ixtest,1), x2val(ixtest,1), x3val(ixtest,1))
        end do

        absinterperror = abs(interpval-actval)
        relinterperror = abs(interpval/actval-1)

        maxabserror = maxval(absinterperror)
        maxrelerror = maxval(relinterperror)

        write (*, *) 'Max absolute interpolation error is', maxabserror
        write (*, *) 'Max relative interpolation error is', maxrelerror
      contains

!The function over which I'll be interpolating
        real (kind=rk) function funcfortestinterp(x1, x2, x3)
          implicit none
          real (kind=rk) :: x1, x2, x3

          funcfortestinterp = x1*x2*(x3)

        end function funcfortestinterp
!---------------------------------------------------------------------------------------------------------!
      end subroutine testlinearinterp3
!---------------------------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------------------------!
subroutine getannuityrates(allyears,survprobs,intrate,annuityrates)

    implicit none

    integer, intent(in) :: allyears
    real (kind=rk), intent(in) :: intrate
    real (kind=rk), intent(in) :: survprobs(allyears)
    real (kind=rk), intent(out) :: annuityrates(allyears)
    
    
    integer :: yearsforpension
    integer :: ixforprob, ixforann
    real (kind=rk), allocatable :: probpay(:)
    
    real (kind=rk) :: recannuityrates(allyears)
    
do ixforann = 1, allyears-1, 1
        yearsforpension = allyears - ixforann 
        
        allocate (probpay(yearsforpension))
        do ixforprob = 1, yearsforpension
       !     probpay(ixforprob) = product(survprobs(ixforann + 1:ixforann+ixforprob))
            probpay(ixforprob) = product(survprobs(ixforann:ixforann+ixforprob - 1))
       
        end do
        
        call getfairvalue(yearsforpension,probpay,intrate,annuityrates(ixforann))
    
        deallocate (probpay)
end do !ixforann   

annuityrates(allyears) = (huge(annuityrates(allyears)))/10.0_rk



end subroutine getannuityrates
! ---------------------------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------------------------!
subroutine getdiscountfrom1(allyears,beta,discountfrom1)

    implicit none

    integer, intent(in) :: allyears

    real (kind=rk), intent(in) :: beta
    real (kind=rk), intent(out) :: discountfrom1(allyears)
    
    integer :: ixyear

    
discountfrom1(1) = 1.0_rk
do ixyear = 2, allyears, 1
   discountfrom1(ixyear) = (beta**(ixyear-1))
end do !ixforprob   

end subroutine getdiscountfrom1
! ---------------------------------------------------------------------------------------------------------!
! ---------------------------------------------------------------------------------------------------------!
subroutine getfairvalue(maxyears,probpay,intrate,annuityrate)
  
  implicit none

  integer, intent(in) :: maxyears
  real (kind=rk) :: probpay(maxyears)
  real (kind=rk) :: intrate
  real (kind=rk) :: annuityrate
  
  !For program
  integer :: ixyear
  real (kind=rk) :: forsum
  
  forsum = 0.0_rk
  do ixyear = 1, maxyears, 1
    forsum = forsum + probpay(ixyear)/((1+intrate)**ixyear) 
  end do
  
  annuityrate = 1.0_rk/forsum
    
end subroutine
! ---------------------------------------------------------------------------------------------------------!


! ---------------------------------------------------------------------------------------------------------!




!---------------------------------------------------------------------------------------------------------!      

! --------------------------------------------------------------------
! --------------------------------------------------------------------
subroutine getmatrixinverse(matrix, dim, dimplus1, inverse)
!Calculate a matrix inverse using the NAG routine. NAG routine gives answer in an odd format - this converts it to a usable format
    use  globalvalues ! only here so that I can access globalrank for an error message
    
   implicit none
        integer, intent(in) :: dim
        integer, intent(in) :: dimplus1
        real (kind=rk), intent(in) :: matrix(dim, dim)
        real (kind=rk), intent(out) :: inverse(dim, dim)
  
        !For program
        real (kind=rk) :: forinverse(dimplus1, dim)
        real (kind=rk) :: copyofinverse(dim, dim)
        integer :: ixi, ixj
        
        !For NAG
        integer :: ifail
        
        !For checking
        real (kind=rk) :: forcheck(dim, dim)
        real (kind=rk) :: sumabs
        
        forinverse = 0.0_rk 
        forinverse(1:dim, :) = matrix
       
       ifail = 1
       call F01ADF(dim, forinverse, dimplus1, ifail)
       if (ifail.eq.0) then
             do ixi = 1, dim
               do ixj = 1, ixi
                 inverse(ixi, ixj) = forinverse(ixi + 1, ixj)
               end do
           end do

          do ixi = 1, dim
               do ixj = ixi, dim
                 inverse(ixi, ixj) = inverse(ixj, ixi)
               end do
          end do
                         
      
          !Check have I got the inverse right
          forcheck = matmul(matrix, inverse) !This should be the identity      
          do ixi = 1, dim
            forcheck(ixi, ixi) = forcheck(ixi, ixi) - 1.0_rk
          end do

          !Get the sum of absolute values
          sumabs = 0.0_rk
          do ixi = 1, dim
              do ixj = 1, dim
                sumabs = sumabs + abs(forcheck(ixi, ixj))
              end do          
          end do

          
          if (sumabs.gt.0.001) then 
                write(*,*) sumabs
              write(*,*) 'error in forinverse - inverse either wrong or innacurate'
              stop
          end if         
        else !        if (ifail.eq.0) then Inverse didn't computer
            inverse = 0 
        write(*,*) 'error in getmatrixinverse, NAG couldnt compute inverse'
        end if
                  
      
  end subroutine getmatrixinverse
  
integer function myMinloc(array, dim) ! The Fortan intrinsic has returned an error
        implicit none
        integer, intent(in) :: dim
        real (kind = rk), intent(in) ::  array(dim)
        
        integer :: ix
        
        myMinloc = 1
        do ix = 2, dim
            if (array(ix).lt.(array(myMinloc))) myMinloc = ix
        end do

end function


! ---------------------------------------------------------------------------------------!
! Extrapolation is by default here. There is no option to use the bottom value (as I have
! in the one-dimensional linear interpolation).

      !linearinterp2_fast has the following difference from linearinterp2
      ! 1. It doesn't call findingr but does it within this module
      ! 2. It checks if the last saved value for the lowerlocs are still correct and if so, doesn't do a binary search
      
      subroutine linearinterp2_fast(x1grid, x2grid, dim1, dim2, x1val, x2val, yval, ygrid)
!        use types
        implicit none
! Arguments
        integer, intent (in) :: dim1
        integer, intent (in) :: dim2
        real (kind=rk), intent (in) :: x1grid(dim1)
        real (kind=rk), intent (in) :: x2grid(dim2)
        real (kind=rk), intent (in) :: x1val
        real (kind=rk), intent (in) :: x2val
        real (kind=rk), intent (out) :: yval
        real (kind=rk), intent (in) :: ygrid(dim1, dim2)
        
        
! For programme
        integer :: x1lowerloc, x2lowerloc
        integer, save :: x1lowerloc_save, x2lowerloc_save
        integer :: scratch_nearestloc
        integer :: newxgridlowerloc
        real (kind=rk) :: newygrid(2)
        real (kind=rk) :: newxgrid(2)
         
        integer ::  x1upperloc
         real (kind=rk) :: xgridmin, xgridmax, xgriddiffinv, weightmax, weightmin
                integer :: jl, ju, jm
       logical ::        binarySearch
integer :: tempCounter, numCycles
        
write(*,*) 'DOESNT SEEM SAFE IN PARALELL'

        !Have commented out for speed
        !if ((size(x1grid).ne.dim1) .or. (size(x2grid).ne.dim2)) then
        !  write (*, *) 'An argument in linearinterp2 is not of the specified length'
        !  stop
        !end if

    !!!___________________________________!!!!
    !!! Find Square
    !!!___________________________________!!!!
                
                
    ! Do find in grid without calling the routine....for speed...calling the routine slows us down
      ! This block of code finds the two points on the xgrid between which the x value lies        
        !call findingridr(x1grid, dim1, x1val, scratch_nearestloc, x1lowerloc)

       binarySearch = .true.
        
       if ((x1lowerloc_save.lt.1).or.(x1lowerloc_save.gt.(dim1 - 1))) then
            binarySearch = .true.
       else if (.not.((x1val.gt.x1grid(x1lowerloc_save)).and.(x1val.lt.x1grid(x1lowerloc_save+1)))) then
            binarySearch = .true.
       else
            binarySearch = .false.
       end if
       
                   
       if (binarySearch) then
                    jl = 0
                    ju = dim1 + 1
                    !do  tempCounter = 1, numCycles1
                    do while (ju-jl.gt.1)
                   
                      jm = (ju+jl)/2
                      if (x1val.gt.x1grid(jm)) then
                        jl = jm
                      else
                        ju = jm
                      end if
                    end do
                   !write(*,*) tempCounter
                    
                    
                    x1lowerloc = jl

                if (x1lowerloc.eq.0) then
                  x1lowerloc = 1
                end if
                if (x1lowerloc.eq.dim1) then
                  x1lowerloc = dim1 - 1
                end if

        else  ! The bounds from last call are still the correct ones      
                
                 x1lowerloc = x1lowerloc_save            
        end if
        

        

        ! Save for next run
        x1lowerloc_save = x1lowerloc

! Find the square
!       call findingridr(x2grid, dim2, x2val, scratch_nearestloc, x2lowerloc)

       binarySearch = .true.
        
       if ((x2lowerloc_save.lt.1).or.(x2lowerloc_save.gt.(dim2 - 1))) then
            binarySearch = .true.
       else if (.not.((x2val.gt.x2grid(x2lowerloc_save)).and.(x2val.lt.x2grid(x2lowerloc_save+1)))) then
            binarySearch = .true.
       else
            binarySearch = .false.
       end if
        
        
      
       if (binarySearch) then

            jl = 0
            ju = dim2 + 1
          
            do while (ju-jl.gt.1)
            !do  tempCounter = 1, numCycles2
              jm = (ju+jl)/2
              if (x2val.gt.x2grid(jm)) then
                jl = jm
              else
                ju = jm
              end if
            end do

            x2lowerloc = jl

            ! Find if any of the values are off the grid
        
              if (x2lowerloc.eq.0) then
                x2lowerloc = 1
              end if
              if (x2lowerloc.eq.dim2) then
                x2lowerloc = dim2 - 1
              end if
        else        
                    ! The bounds from last call are still the correct ones
                    x2lowerloc = x2lowerloc_save
        end if        
        
       

          ! Save for next run
          x2lowerloc_save = x2lowerloc


 !!! Do this but faster         
        !call linearinterpfromlocations(dim1, x1grid, ygrid(:,x2lowerloc), x1lowerloc, x1val,  1, 1, newygrid(1))
         
       if (x1lowerloc.ne.0 .and. x1lowerloc.ne.dim1) then
          x1upperloc = x1lowerloc + 1
          xgridmin = x1grid(x1lowerloc)
          xgridmax = x1grid(x1upperloc)
          xgriddiffinv = 1/(xgridmax-xgridmin)
          weightmax = (x1val-xgridmin)*xgriddiffinv
          weightmin = (xgridmax-x1val)*xgriddiffinv
          
          ! Calculating the interpolated function values
          newygrid(1) = weightmin*ygrid(x1lowerloc, x2lowerloc) + weightmax*ygrid(x1upperloc,x2lowerloc)
       else
           call linearinterpfromlocations_fast(dim1, x1grid, ygrid(:,x2lowerloc), x1lowerloc, x1val,  1, 1, newygrid(1))
       end if      
          
          
          
  !!! Do this but faster        
!          call linearinterpfromlocations(dim1, x1grid,ygrid(:,x2lowerloc+1), x1lowerloc, x1val,  1, 1, newygrid(2))

       if (x1lowerloc.ne.0 .and. x1lowerloc.ne.dim1) then
          x1upperloc = x1lowerloc + 1
          xgridmin = x1grid(x1lowerloc)
          xgridmax = x1grid(x1upperloc)
          xgriddiffinv = 1/(xgridmax-xgridmin)
          weightmax = (x1val-xgridmin)*xgriddiffinv
          weightmin = (xgridmax-x1val)*xgriddiffinv
          
          ! Calculating the interpolated function values
          newygrid(2) = weightmin*ygrid(x1lowerloc, x2lowerloc + 1) + weightmax*ygrid(x1upperloc,x2lowerloc + 1)
       else
           call linearinterpfromlocations_fast(dim1, x1grid, ygrid(:,x2lowerloc), x1lowerloc, x1val,  1, 1, newygrid(2))
       end if      
          
          
          
          
! Now generate a new grid from from the x2 bounds (a new ygrid has already bee generated)
        newxgrid(1) = x2grid(x2lowerloc)
        newxgrid(2) = x2grid(x2lowerloc+1)

!Think this is redundant
        if (x2lowerloc.eq.0) then
          newxgridlowerloc = 0
        else if (x2lowerloc.eq.dim2) then
          newxgridlowerloc = 2
        else
          newxgridlowerloc = 1
        end if
        call linearinterpfromlocations_fast(2, newxgrid, newygrid, newxgridlowerloc, x2val,  1, 1, yval)

      end subroutine linearinterp2_fast

      !linearinterpfromlocations_fast just has reordered arguments (why I do not know - I did this in optimality) 
      subroutine linearinterpfromlocations_fast(length, xgrid, ygrid, xlowerloc, xval, extrapbot, extraptop, yval)
        implicit none
        integer, intent (in) :: length
        real (kind=rk), intent (in) :: xgrid(length)
        real (kind=rk), intent (in) :: ygrid(length)
        integer, intent (in) :: xlowerloc
        real (kind=rk), intent (in) :: xval
        integer, intent (in) :: extrapbot
        integer, intent (in) :: extraptop
        real (kind=rk), intent (out) :: yval


! For programme
        real (kind=rk) xgridmin
        real (kind=rk) xgridmax
        real (kind=rk) xgriddiffinv
        real (kind=rk) weightmin
        real (kind=rk) weightmax
        real (kind=rk) extrapslope
        integer :: xupperloc

! First, and easiest, get the interpolated value if xval is interior to the grid
        if (xlowerloc.ne.0 .and. xlowerloc.ne.length) then
          xupperloc = xlowerloc + 1
          xgridmin = xgrid(xlowerloc)
          xgridmax = xgrid(xupperloc)
          xgriddiffinv = 1/(xgridmax-xgridmin)

          weightmax = (xval-xgridmin)*xgriddiffinv
          weightmin = (xgridmax-xval)*xgriddiffinv

! Calculating the interpolated function values
          yval = weightmin*ygrid(xlowerloc) + weightmax*ygrid(xupperloc)
        else if (xlowerloc.eq.0) then
          if (extrapbot.eq.0) then
            yval = ygrid(1)
          else
            extrapslope = (ygrid(2)-ygrid(1))/(xgrid(2)-xgrid(1))
            yval = ygrid(1) + extrapslope*(xval-xgrid(1))
          end if

        else ! a is hu-eg so we extrapolate
          if (extraptop.eq.0) then
            yval = ygrid(length)
          else
            extrapslope = (ygrid(length)-ygrid(length-1))/(xgrid(length)-xgrid(length-1))
            yval = ygrid(length) + extrapslope*(xval-xgrid(length))
          end if
        end if


      end subroutine linearinterpfromlocations_fast

      
      
            subroutine linearinterp1_fast(xgrid, ygrid, length, xval, yval, extrapbot, extraptop)
! Subroutine for linear interpolation
! Arguments are 1. (input) xvlalues
! 2. (input) yvalues
! 3. (input) length of xgrid and ygrid,
! 4. (input) value of x to interpolate at
! 5. (output) interpolated value
! 6. (input) indicator for whether to extrapolate below bottom point of grid
! 7. (input indicator for whether to extrapolate above top point of grid

! If the choose not to extrapolate interpolated value is set equal to either the top or the bottom
! point of the grid
! It is fully generic

        implicit none
! Arguments
        integer, intent (in) :: length
        integer, intent (in) :: extrapbot
        integer, intent (in) :: extraptop
        real (kind=rk), intent (in) :: xgrid(length)
        real (kind=rk), intent (in) :: ygrid(length)
        real (kind=rk), intent (in) :: xval
        real (kind=rk), intent (out) :: yval

! For programme
        integer :: xlowerloc
        integer, save :: xlowerloc_save
        integer scratch_xnearestloc ! This areguments required by findingridr but I'm not using them here - hence
        real (kind=rk) xgridmin
        real (kind=rk) xgridmax
        real (kind=rk) xgriddiffinv
        real (kind=rk) weightmin
        real (kind=rk) weightmax
        real (kind=rk) extrapslope
        integer :: xupperloc
        logical :: binarySearch

        
        ! scratch
        integer :: jl, ju, jm
        intrinsic size

        
! Checking inputs - have commented this out for speed
     !   if ((size(xgrid).ne.length) .or. (size(ygrid).ne.length)) then
     !     write (*, *) 'An argument in linearinterp1 is not of the specified length'
     !     stop
     !   end if

     !   if ((extrapbot.ne.0 .and. extrapbot.ne.1) .or. (extraptop.ne.0 .and. extraptop.ne.1)) then
     !     write (*, *) 'An argument in linearterp1 that needs to be either 0 or 1 is not so'
     !     stop
     !   end if

        if (abs(xgrid(length)-xgrid(1)).lt.0.1d-15) then ! If the top of the grid is equal to the bottom
          yval = ygrid(1) ! (which would include the scenario that all the entries are the same) then I
! assign the bottom value of y to yval
        else
      
       binarySearch = .true.
        
       if ((xlowerloc_save.lt.1).or.(xlowerloc_save.gt.(length - 1))) then
            binarySearch = .true.
       else if (.not.((xval.gt.xgrid(xlowerloc_save)).and.(xval.lt.xgrid(xlowerloc_save+1)))) then
            binarySearch = .true.
       else
            binarySearch = .false.
       end if
             
            
        if (binarysearch) then        
        ! Find the bounds    
        jl = 0
        ju = length + 1
        do while (ju-jl.gt.1)
          jm = (ju+jl)/2
          if (xval.gt.xgrid(jm)) then
            jl = jm
          else
            ju = jm
          end if
        end do

        xlowerloc = jl
            
        else  ! The bounds from last call are still the correct ones      
                
                 xlowerloc = xlowerloc_save            
        end if
        

        

        ! Save for next run
        xlowerloc_save = xlowerloc
    
        
          !call findingridr(xgrid, length, xval, scratch_xnearestloc, xlowerloc)
          !call linearinterpfromlocations_fast(length, xgrid, ygrid, xlowerloc, xval,  extrapbot, extraptop, yval)
         
! First, and easiest, get the interpolated value if xval is interior to the grid
        if (xlowerloc.ne.0 .and. xlowerloc.ne.length) then
          xupperloc = xlowerloc + 1
          xgridmin = xgrid(xlowerloc)
          xgridmax = xgrid(xupperloc)
          xgriddiffinv = 1/(xgridmax-xgridmin)

          weightmax = (xval-xgridmin)*xgriddiffinv
          weightmin = (xgridmax-xval)*xgriddiffinv

! Calculating the interpolated function values
          yval = weightmin*ygrid(xlowerloc) + weightmax*ygrid(xupperloc)
        else if (xlowerloc.eq.0) then
          if (extrapbot.eq.0) then
            yval = ygrid(1)
          else
            extrapslope = (ygrid(2)-ygrid(1))/(xgrid(2)-xgrid(1))
            yval = ygrid(1) + extrapslope*(xval-xgrid(1))
          end if

        else ! a is hu-eg so we extrapolate
          if (extraptop.eq.0) then
            yval = ygrid(length)
          else
            extrapslope = (ygrid(length)-ygrid(length-1))/(xgrid(length)-xgrid(length-1))
            yval = ygrid(length) + extrapslope*(xval-xgrid(length))
          end if
        end if
        
        end if

            end subroutine linearinterp1_fast

   logical function readinlogical(realoneorzero, whichvalinallpoints, oneistrueallothersfalse)
          implicit none
          
          real (kind=rk) :: realoneorzero
          integer :: whichvalinallpoints
          logical, optional :: oneistrueallothersfalse
          logical :: inprogoneistrueallothersfalse
          
          if (present(oneistrueallothersfalse)) then
          	inprogoneistrueallothersfalse = oneistrueallothersfalse
          else 
          	inprogoneistrueallothersfalse = .false.
          end if
          
          if (inprogoneistrueallothersfalse) then
             if ((realoneorzero.gt.0.999_rk) .and. (realoneorzero.lt.1.001_rk)) then
                readinlogical  = .true.
              else
                readinlogical = .false.
              end if
          else     
              if ((realoneorzero.gt.-0.001_rk) .and. (realoneorzero.lt.0.001_rk)) then
                readinlogical = .false.
              else if ((realoneorzero.gt.0.999_rk) .and. (realoneorzero.lt.1.001_rk)) then
                readinlogical  = .true.
              else
                write (*, *) 'error in readindata - an element should be either zero or 1'
                write (*, *) 'element is', whichvalinallpoints
                stop
              end if
          end if
          
   end function readinlogical
            
   
   
    end module routines_generic
